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SUMMARY 


This paper discusses a research effort directed toward increased efficiency in 
calculating transient temperature fields in complex aerospace vehicle structures. 
Explicit solution techniques which require minimal computation per time step and 
implicit techniques which permit larger time steps because of better stability are 
reviewed and evaluated. A set of implicit solution algorithms with variable time 
steps (GEARIB) is described and evaluated. Test problems for evaluating the algo- 
rithms include a coarse model of the Space Shuttle Orbiter wing, an insulated frame 
test article, a metallic panel for a thermal protection system, and detailed models 
of single-bay and three-bay sections of the Shuttle wing. Results generally indicate 
that implicit algorithms, especially GEARIB with variable step size and order, are 
faster than explicit algorithms for transient structural heat-transfer problems when 
the governing equations are stiff . Stiff equations occur in many practical problems 
particularly in insulated thin metal structures. 

The effect on algorithm performance of different models of an insulated cylinder 
was also studied. The study revealed that the stiffness of the problem is highly 
sensitive to modeling details and that careful modeling can reduce the stiffness of 
the resulting equations to the extent that explicit methods may be advantageous. 

Studies of partitioning techniques for improving the performance of solution 
algorithms were carried out. A mixed implicit-explicit technique developed for 
dynamic response has been adapted to thermal problems and is demonstrated. In this 
method, the model is separated into stiff and nonstiff regions. The implicit algo- 
rithm is applied to the stiff portion and the explicit to the nonstiff portion. Two 
operator- splitting techniques for speeding up the solution of the algebraic equations 
associated with implicit algorithms are also described and demonstrated. Both are 
based on separating the coefficient matrices into two parts and solving the resulting 
equations by iteration. The first technique bases the separation on stiff and non- 
stiff elements of the structure, and the second bases the separation on the sparsity 
structure of the matrix. The most effective of the techniques tried was a separation 
based on matrix sparsity (using an incomplete Cholesky decomposition) combined with 
conjugate gradient iteration. Especially noteworthy is the fact that the performance 
of this technique is insensitive to the band structure of the matrix. 


INTRODUCTION 

An effort is in progress at the Langley Research Center to gain increased effi- 
ciency in the prediction and optimization of the thermal-structural behavior of aero- 
space vehicle structures. A principal task is to reduce the computing effort for 
obtaining transient temperatures. This paper is focused on (1) assessment of the 
performance of explicit and implicit temporal integration algorithms, (2) the effects 
of modeling on the performance of the algorithms, and (3) techniques for solving the 
algebraic equations associated with implicit methods. 

In the current and recent literature, many of the difficulties associated with 
the solution of transient heat- transfer problems and other time- dependent physical 



problems are associated with the stiff ^ nature of the governing differential equa- 
tions (refs. 1 to 7). Stiff sets of equations for thermal applications occur in 
problems where the rate of heating is significantly slower than the speed of 
propagation of temperature differences between adjacent points in the model. In 
particular, insulated thin metal structures under long-duration heating often lead to 
stiff equations. A preference is evident among researchers for implicit algorithms 
for solution of stiff sets of ordinary differential equations. However, many 
engineering analysts prefer to use the longer established explicit algorithms, even 
for stiff problems. A partial explanation for this dichotomy is that the full power 
of the implicit approach has not been transferred from researchers to engineering 
analysts. In particular, implicit algorithms are usually implemented in computer 
programs with a fixed time step (refs. 8 to 10). 

In the explicit algorithms, the time step is limited (often severely) in order 
for the technique to be stable. In the implicit algorithms, there is no stability- 
imposed limitation on step size, and the step size is limited by solution accuracy 
only. This conclusion holds in general for linear systems, but it is also found to 
hold in many nonlinear applications (ref. 11). Thus implicit algorithms can, in 
general, use much larger time steps than explicit algorithms; this is especially true 
for stiff problems. Because a single explicit time step is computationally faster 
than a single implicit time step, the key to the advantageous use of implicit 
algorithms is to use the largest possible time step size. 

The strategy advocated for the solution of problems by implicit methods is to 
use algorithms with variable step size and order and to automatically select both 
throughout the solution process (refs. 12 to 15). A promising set of algorithms, 
developed to implement this strategy, is denoted the GEAR algorithms (refs. 13 and 
14). A version of the GEAR algorithms well-suited to heat-transfer analysis, denoted 
GEARIB, has been recently installed in the SPAR finite-element thermal analyzer 
(ref. 8) for testing. Early evaluations of GEARIB in SPAR and comparisons with other 
implicit algorithms were described in reference 16 for the first three problems in 
this paper. 

The first objective of the present paper is to describe recent evaluations, 
improvements, and demonstrations in the use of explicit and implicit algorithms for 
•transient thermal analysis of structures. A coarse model of the Space Shuttle 
Orbiter wing, an insulated frame test article, a metallic multiwall thermal protec- 
tion system panel, and detailed models of single-bay and three-bay sections of the 
Shuttle wing are analyzed. Comparisons between implicit and explicit algorithms are 
presented and the performance of the GEARIB algorithms, especially the value of vari- 
able step size, is demonstrated. 

A characteristic of thermal analysis by finite-element and lumped-parameter 
techniques is that modeling strongly affects the stiffness. Since stiffness is one 
of the key factors in the performance of implicit and explicit algorithms, the second 
objective of the paper is to study the effects of modeling. This paper describes a 
study carried out for an insulated cylinder of the effects of modeling on the perfor- 
mance of explicit and implicit algorithms. 
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Stiff sets of ordinary differential equations are characterized by solutions 
with widely varying time constants. The typical case is when the solution to the 
homogeneous problem has some very small time constants compared with those of the 
forcing function (ref. 1). 
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Finally, when the stiffness of the problem is due to only part of the finite- 
element model, partitioning techniques may be useful. The third objective of the 
paper is to explore the potential of some partitioning techniques. The first is a 
mixed implicit- explicit algorithm developed by Hughes and Liu (refs. 17 and 18) which 
was adapted to thermal analysis in reference 19 and is demonstrated for an insulated 
panel. Also, two techniques for speeding up the solution of the algebraic equations 
associated with implicit algorithms are described and demonstrated. Both are based 
on separating the coefficient matrices into two parts and solving the resulting equa- 
tions by iteration. The first technique bases the separation on stiff and nonstiff 
elements of the structure and the second bases the separation on the sparsity 
structure of the matrix (ref. 20). 


SYMBOLS 

C capacitance matrix 

c heat capacity 

Jr 

DT integration time step size in SPAR program (ref. 8) 

F right-hand side of equations for transient temperature 

h n nth time step 

K conductivity matrix 

k thermal conductivity 

Q thermal load vector 

R residual of system of equations generated by implicit method 

T vector of temperatures; dot over T indicates differentiation with 

respect to time 

t time 

tg typical time constant of applied heat loads 

t nth time point 

t s stability limit of time step for explicit algorithms 

a thermal diffusivity 

coefficients in backward-difference method (eq. (6)) 

coefficients in backward-difference and Adams-Moulton methods (eqs. (5) 
and ( 6 ) ) 

p mass density 
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Subscript; 


n time step number 


NATURE OF ALGORITHMS USED IN TRANSIENT STRUCTURAL THERMAL ANALYSIS 

The governing semi-discrete equations for a transient heat-transfer problem 
discretized by a finite-element or finite-difference technique are as follows: 


CT = Q(T,t) - K(T,t) T = F(T/t) (T(0) - T ) 


( 1 ) 


where T is a vector of temperatures, C and K are matrices, Q and F are 
vectors, and F is generally a nonlinear function of T. Obtaining an analytical 
solution to equation ( 1 ) is usually impractical and numerical integration methods are 
used. The simplest numerical integration technique is the Euler method which uses 
the first two terms in a Taylor series to predict T 


T(t n + 1> = T( V 


+ h T(t ) = T(t ) + h C -1 f[" T( t ),t 1 
n n n n L n nj 


(2) 


where h n is the time step size at the nth step. Euler's method is an example 
of an explicit integration technique, so- named because T(t n+1 ) is given explicitly 
in terms of known quantities. Another approach to the numerical integration of 
equation (1) is the backward-dif f erence method which is an example of an implicit 
method. In this approach. 


T(t n + 1 > - T( V + h „ *<W - T( V + v ' 1 F [ T( W'Vi] (3 > 

Equation (3) is a system of implicit equations for T(t n+1 ), which is generally non- 
linear. The explicit algorithm is therefore easier to implement; but, the time step 
must be bounded to avoid numerical instability (unbounded propagation of numerical 
errors during the solution) . Implicit techniques are generally unconditionally 
stable for linear systems (and some nonlinear applications) and thus can take larger 
time steps which are determined from accuracy considerations rather than stability. 

For one- dimensional conduction problems, the time for a temperature disturbance 
to propagate from one point to another is roughly proportional to L 2 /a, where L is 
the distance between the two points and a is the thermal diffusivity. The time 
step t required to insure stability of an explicit integration method is propor- 
tional to this intrinsic time. For a finite-element discretization, the explicit 
time step is determined by the element with the smallest value of L 2 /a, where L is 
the length of the element. The rapidity of variation of the temperature, on the other 
hand, depends on the time scale of the applied heating. Accuracy considerations 
require a time step which is comparable to a typical time constant of the applied 
heat loads tg. By definition, a problem is stiff whenever tg is much larger 
than t g . Implicit methods (in which time steps are bounded only by accuracy 
requirements) can use large time steps (of order tg) , whereas explicit methods 
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which are bound by stability requirements require small time steps (of order 
L^/oc). It follows that stiff problems are usually best solved by implicit methods. 
The effort involved in solving a system such as equation (3) is usually cost- 
effective if a small number of large time steps are used. The preceding comments 
hold for general transient thermal problems except that the intrinsic time steps have 
more complex forms than L^/a. 


The Euler method and the backward-difference methods are presented as examples 
from a large class of explicit and implicit techniques, respectively. Higher order 
methods (e.g. , multistep) in both the explicit and implicit classes typically use 
more previous information to predict the temperature but the stability properties 
of explicit multistep methods are similar to those of the Euler method (ref. 11). 
Explicit methods are unstable for time steps much larger than the smallest value of 
L ^/a. Accordingly, thermal analysis computer programs generally select explicit 
time steps automatically based on the stability requirement. When using implicit 
methods, the analyst is left to select the implicit time step and order without a 
great deal of guideline information, and often several trial runs are needed. There 
is an emerging consensus that the preferred approach for integrating stiff systems 
of ordinary differential equations is to use implicit methods which automatically 
select the order and the step size based on desired accuracy. One software package, 
denoted the GEARIB algorithms, has these features and is discussed subsequently. 


THE GEARIB ALGORITHMS 

Several software packages based on the work of Gear have been developed for 
general use (ref. 13). The package most appropriate for application to finite- 
element thermal analysis is denoted GEARIB. This package is intended to solve 
systems of ordinary differential equations of the form 


CT = F(T,t) 


(4) 


The package employs two classes of implicit multistep methods, Adams-Moulton and 
general backward differences. For nonstiff equations the Adams-Moulton method of 
orders 1 through 12 is used. This method has the general form 

q 

T(t = T(t ) + h > p. T(t _ .) (5) 

n+1 n n / j i n+1-i ' 

i=0 

where q is the order. For stiff equations, backward-difference algorithms of 
orders 1 through 5 are used. These have the general form 

q 

T(t n + 1> - Vo f <W + Y, “i T(t n + ,-i> <«> 

i=1 

The coefficients a ^ and (3^ are given in reference 15. GEARIB employs a 
predictor-corrector approach where equation (5) or (6) is the corrector equation and 
extrapolation is used for the predictor. The user selects the class of methods 
(Adams-Moulton or backward differences) and error tolerance. GEARIB automatically 
selects the step size and order by a technique described in reference 13. The 
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essence of the step size and order selection technique is that at each time step the 
error tolerance is compared with an estimate of the actual error resulting from use 
of the previous step size and order. Both are scaled to force the actual error to be 
no larger than the specified error tolerance. 

Use of the GEARIB algorithms is illustrated by using the backward- difference 
option (eq. (6)). Applied to equation (4), equation (6) gives 


R = C 


T(t „ + ,> 


-E 


i=1 


o. T(t . ) 
l n+1-i 


- h B n 
n 0 


T(t „ + 1>' t n + 1 


■]- 


(7) 


This system of nonlinear algebraic equations is solved by the modified Newton 
method; that is, 

Ti+1<t n+1 > ■ Ti «Vi> - [i] ’* 

where 


5R 

ar 


C 


flh J 
0 n 


and J = SF/9T is the Jacobian of the system. Methods and options used in GEARIB 
for computing J are described in references 13 and 21. 


DESCRIPTION OF TEST PROBLEMS AND RESULTS 

As discussed previously, the stiffness of a transient thermal analysis problem 
is in part dependent on the size of the finite elements used in the modeling. For 
the same heat loads, a detailed model which has small finite elements will be stiff er 
than a coarse model which has larger elements. In assessing the performance of 
explicit and implicit algorithms, it is important to include examples which span the 
range of model fineness used in practical applications. The test problems were 
chosen with this in mind. The Space Shuttle Qrbiter wing model represents a coarse 
model suitable for overall trends in early Shuttle thermal analyses. The multiwall 
thermal protection system panel represents, probably, the most refined model that 
would be used. The other three examples represent typical models which are somewhere 
in between the aforementioned two extremes. 

In comparing the efficiency of the different algorithms, the Euler method was 
used as the only explicit method, whereas the first-order backward-difference, the 
Crank-Nicholson, and the GEARIB algorithms were used to represent implicit methods. 
Euler was chosen because, of all the power series explicit algorithms available in 
SPAR, the Euler algorithm (based on the first two terms of the power series) is the 
most efficient for stability-limited problems. The stability limit for the time step 
in SPAR is estimated as t g = min where c^ and k^ are the diagonal 

terms of the capacitance and conductance matrices, respectively. The conductance 
matrix includes radiation effects. The actual time step used in the analysis was 
selected to be 95 percent of the estimated stability limit, unless otherwise noted. 
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Solution times reported in comparing the algorithms are central processor unit (CPU) 
times. (A comparison of solution times for all five examples is presented in 
table 1.) All calculations were performed on the Control Data CYBER 173 computer 
system at the Langley Research Center under network operating systems (NOS) 1.3 and 
1.4 unless otherwise noted. It is recognized that the use of GEARIB to represent 
implicit algorithms and the Euler method to represent explicit algorithms is somewhat 
biased. However, for stability-limited parabolic problems such as the stiff thermal 
analysis problems herein, the use of more advanced explicit algorithms are not 
expected to be significantly more efficient than the Euler method. Thus, the bias is 
probably not very significant. 


Space Shuttle Orbiter Wing 

The Space Shuttle Orbiter wing (fig. 1(a)) was chosen as the first example 
problem. The SPAR finite-element model (fig. 1(b)) consists of a relatively coarse 
idealization of the structure augmented by a representation of the insulation 
attached to the upper and lower surfaces. This model is useful for obtaining rough 
overall qualitative ideas about the temperature distributions but is too coarse for 
a detailed thermal analysis. The structure is modeled by one-dimensional, triangu- 
lar, and quadrilateral (K21, K31, and K41, respectively) SPAR conduction elements. 

The insulation on each surface is modeled by six layers of one-dimensional conduction 
elements. Use of these latter elements neglects lateral heat transfer in the insula- 
tion. This modeling decision was guided by results from a previous wing model 
(ref. 21) in which lateral heat transfer in the insulation was included. It was 
found that lateral temperature gradients were negligible when compared with gradients 
through the insulation. The complete model contains 2289 grid points, 840 one- 
dimensional ( 1-D) and 560 two-dimensional (2-D) elements in the structure, and 1962 
one- dimensional elements in the insulation. Thermal properties of the aluminum 
structure are temperature dependent ; thermal properties of the insulation are temper- 
ature and pressure dependent. The pressure dependence is treated in SPAR as time 
dependence using the variation of pressure as a function of time from the trajectory 
data for a simulated flight. As a result of these dependencies, the conductivity of 
the insulation in the outer layers changes by a factor of nine during the temperature 
history. 

For the purpose of this analysis, the heating on the wing is represented by a 
specified temperature history on the external surface of the insulation on the under 
side of the wing (shown as the solid line in fig. 2) and is roughly indicative of 
atmospheric reentry heating. Temperature histories of the wing for 4500 seconds were 
computed by using the explicit (Euler), implicit (Crank-Nicholson) , and GEARIB algo- 
rithms. Figure 2 shows the temperature histories of a point in the structure and a 
point in the insulation at a typical wing cross section. The explicit, implicit ,' and 
GEARIB algorithms produced essentially the same temperature histories. Solution time 
comparisons are shown in table 1(a) along with the time steps used to obtain compa- 
rable accuracy. The explicit algorithm used a time step of 10 seconds - in fact, the 
stability limit calculated by SPAR (t g ) was over 100 seconds but the time step size 
was dictated by accuracy and the need to periodically update temperature-dependent 
material properties. The large permitted time step is due to the coarse modeling of 
the structure which did not include any thin, high- conducting, or radiating elements 
which can lead to stiff equations. The implicit algorithm (Crank-Nicholson) required 
a time step of 10 seconds to achieve comparable accuracy and required about five 
times as much computer time as the ejq>licit algorithm. The GEARIB algorithms 
performed very well for this problem. By adaptively varying the time step from 
1.0 second early in the temperature history to as large as 528 seconds toward the 
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end, GEARIB required only 570 CPU seconds to complete the solution. GEARIB required 
only eight recalculations and refactorings of the Jacobian, which shows that even 
with the large variations in material properties, the problem is only mildly 
nonlinear. 


Insulated Frame Test Article 

An insulated frame test article analyzed and tested under transient heating as 
described in reference 22 is shown in figure 3(a). The test article consisted of an 
aluminum skin/stringer structure with two corrugated aluminum frames. An installa- 
tion of FRSI (flexible reusable surface insulation) was placed on the skin surface. 
The test article was also equipped with an auxiliary insulation blanket to prevent 
direct heating of the test article sides and back. The principal purpose of the 
study of the configuration, as discussed in reference 22, was to evaluate the thermal 
performance of Shuttle FRSI during a simulated flight. A secondary purpose was to 
evaluate by comparison with test data the adequacy of thermal analysis techniques 
needed by the Shuttle contractors for preflight thermal and thermal stress analysis. 

Because of symmetry, the lumped- parameter model from reference 22 (not shown 
in the figure) consisted of a two-dimensional section of half the structure. The 
lumped-parameter model was converted to a finite- element model for analysis with 
SPAR. The SPAR finite-element model (fig. 3(b)) contains 149 grid points and 
148 elements including one- dimensional elements which account for conduction in the 
aluminum structure and radiation across the air gap and two-dimensional elements 
which model conduction in the insulation and across the gap. Minor modifications 
were made to the finite-element model following the conversion from the lumped- 
parameter model. These modifications consisted of eliminating or consolidating some 
extremely thin or short finite elements in the aluminum structure in order to reduce 
the stiffness of the equations and to increase the allowable time step for the 
ejq>licit solution algorithm. The properties of the aluminum structure are functions 
of temperature, and, as in the previous example, the properties of the insulation are 
functions of temperature and pressure. The heating is simulated by a temperature 
history at the outer surface of the insulation. 

The temperature history for the frame was computed by using the explicit (Euler) 
and implicit techniques ( Crank-Nicholson and backward differences) and GEARIB. Com- 
parisons of solution times are given in table 1(b) . The explicit algorithm had a 
stability limit t of 0.16 second and required 1723 seconds. This time step was 
controlled by conduction through the aluminum elements along the center and front of 
the model. The solution times using the Crank-Nicholson algorithm varied from 475 to 
65 seconds corresponding to variations in the time step from 1 . 0 and 50 seconds . The 
solution times for backward differences were essentially the same as those of Crank- 
Nicholson and are not shown. The GEARIB algorithm used time steps from 50 to 
170 seconds and the solution time was 54 seconds. As indicated in table 2, there is 
very little loss of accuracy in either the structure or insulation temperatures with 
increased time- step size. 

The accuracy of the solutions obtained by the various techniques is further 
assessed in figure 4 which displays temperature histories at a point in the outer 
layer of the aluminum structure corresponding to node 309. (See fig. 3.) The solid 
curve in figure 4 represents specified temperatures at the outer surface of the insu- 
lation (node 29) . The dashed line shows temperatures obtained by the SPAR analyses 
which are plotted as a single curve since there is little difference between the 
results. The long-dash — -short-dash curve shows analytical results from the lumped- 
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parameter analysis of reference 22, which are also in close agreement with the SPAR 
temperatures. The circles in figure 4 represent test data from reference 22. The 
close agreement between the numerical and experimental results indicates that the 
models are adequate to simulate the temperature history in the test article. 


Multiwall Thermal Protection System Panel 

The next example problem is one which evolved out of a study of the thermal 
performance of a titanium multiwall thermal protection system (TPS) panel which is 
under study for future use on space transportation systems (ref. 23). A high level of 
modeling detail was required to accurately represent the local thermal conductivity 
of the multiwall concept and to assess some approximations used in the work reported 
in reference 23. The configuration, depicted in figure 5(a), consists of alternating 
layers of flat and dimpled sheets fused at the crests of the dimpled sheets to form a 
sandwich. The representation of a typical dimpled sheet is shown in figure 5(b) . 

For the purpose of this analysis, it is assumed that the heating does not vary later- 
ally. This assumption, in addition to the regular geometry of the structure, leads 
to the modeling simplification wherein only a triangular prismatic section (indicated 
in fig. 5(a)) of the panel needs to be modeled. The intersection of this prism with 
a typical dimpled layer is indicated by the shaded triangle in figure 5(b). 

The finite-element model shown in figure 5(c) contains 333 grid points located 
on nine titanium sheets (five horizontal and four inclined). The model contains 288 
triangular and quadrilateral metal conduction elements, 264 solid elements (K61 and 
K81) which account for air conduction between the layers, and 544 triangular and 
quadrilateral radiation elements which account for radiation heat transfer between 
adjacent horizontal and inclined sheets. Thermal properties of titanium and air are 
functions of temperature. Radiation exchange (view) factors were computed and sup- 
plied to SPAR with the TRASYS II computer program (ref. 24). 

The temperature history in the panel resulting from a specified transient tem- 
perature at the outer surface of the panel was computed for 3200 seconds. Results 
were obtained with SPAR by using explicit, Crank-Nicholson, backward-difference, and 
GEARIB algorithms. Solution time comparisons are presented in table 1(c). The 
esqplicit algorithm had a stability- limited time step t_ of 0.007 second. This time 
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step was dictated by the short conduction paths between the vertices of adjacent 
triangular layers. The small time step indicate.s that this is an extremely stiff 
problem as may be expected because of the fine detail in modeling. Required solution 
time for the explicit algorithm was estimated by extrapolation to be 98 400 seconds. 

The Crank-Nicholson solution was carried out with time steps of 1.0 and 
5 seconds, which led to solution times of 28 400 (estimation based on extrapolation) 
and 6352 seconds, respectively. Backward difference was used with the same time 
steps and had the same solution times . GEARIB took time steps ranging between 1 . 0 
and 113 seconds and required a solution time of 2754 seconds. A plot of typical 
temperature histories for a point midway through the panel and the primary structure 
are shown in figure 6 along with the applied outer surface temperature. The results 
shown are from the implicit algorithm with a time step of 5 seconds and are identical 
to results with a time step of 1 second and GEARIB. 
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Three Bays of Shuttle Wing 


A three-bay section of the Shuttle wing (fig. 7) was analyzed by the authors 
of reference 25 and the model was made available for the present study. The finite- 
element model shown in figure 8 for wing station 328, unlike the coarse model of the 
complete wing shown previously, is typical of the level of detail used in Shuttle 
thermal modeling practice. The purpose of the study described in reference 25 
was to use thermal analysis tools to predict temperatures measured during the first 
Shuttle flight. Temperatures predicted by the SPAR program were found to be in good 
agreement with flight data. The forward spar web of bay 1 is made of aluminum honey- 
comb core sandwich plate, and the rest are corrugated plates. The upper and lower 
skins are hat-stringer stiffened. The lower skin is covered with HRSI (high- 
temperature reusable surface insulation). The upper skin is covered with LRSI (low- 
temperature reusable surface insulation) and FRSI. 

In order to account for the spanwise heat flow and the effect of the rib 
trusses, the wing segment was modeled in three dimensions. The finite-element model 
has 915 joint locations. The modeling was limited to the major load-carrying 
portion. The wing skins and spar webs, rib caps, shear webs, insulation surface 
coatings, and RTV (room temperature vulcanized) adhesive layers lying on both sides 
of the SIP (strain isolator pad) were modeled with SPAR K41 elements. The spar caps, 
rib caps, and rib trusses were modeled with SPAR K21 elements. The insulation was 
modeled in 13 layers on the lower surface and 5 layers on the upper surface with SPAR 
K81 elements. The SIP was modeled by one layer of SPAR K81 elements. Aerody- 
namically heated surfaces were modeled with one layer of K41 elements. For the 
external and internal radiant heat energy exchanges, a layer of SPAR R41 elements 
(four-node radiation exchange elements) was attached to the outer surface of the 
insulation and the exposed aluminum surfaces. For the spar webs with two sides 
exposed (lying between the bays), one layer of R41 elements was used for each exposed 
surface. Outer surfaces of the forward and rear spar webs were totally insulated. 
Internal convection and external convective cooling (negligible during entry) were 
ignored. The heating rates vary along the model from bay to bay. The heating rates 
used are from reference 25 and are shown in figure 9. 

A 3500- second history of temperature for the model was computed by use of the 

explicit (Euler), implicit (Crank-Nicholson) , and GEARIB methods. Solution time 

comparisons are given in table 1(d). The rather nonstiff nature of this model leads 

to a relatively large explicit time step t of 2.6 seconds. The Crank-Nicholson 

s 

algorithm required a time step of 5 seconds to obtain comparable accuracy. Not 
surprisingly, the explicit solution time was significantly less than the Crank- 
Nicholson solution time. The GEARIB method used time steps varying from 0.1 to 
229 seconds for a solution time that is less than one-fifth the time required by 
the explicit technique and it obtained equal accuracy. 


Single Bay of Shuttle Wing 

The single-bay model shown in figure 10 represents a two-dimensional section 
through the wing in the third bay at wing station 240. The model was obtained from 
the authors of reference 25. The SIP and RTV layers (lying on both sides of the SIP) 
were modeled with SPAR K41 elements. The spar webs and caps, rib caps, and trusses 
and skins were modeled with SPAR K21 elements so as to form a frame. The insulation 
was modeled in 10 layers on the lower surface and 3 layers on the upper surface with 
SPAR K41 elements. Radiation heat transfer from the upper and lower surfaces and 
inside the cavity was modeled with R21 elements. The single-bay model contained a 
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total of 123 grid points and 191 elements. The heat rates as functions of time were 
obtained from the authors of reference 25 and are shown in figure 11. 

A 3500-second temperature history of the model was performed by use of the 
Euler, Crank-Nicholson, and GEARIB algorithms. Solution time comparisons are given 
in table 1(e). The model is stiff due to the short heat paths through the aluminum 
elements near the corner of the frame. The time step needed by the Euler technique 
t„ was 0.1 second. The Crank-Nicholson method was used with a time step of 
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1 .0 second resulting in a solution time which was one-third that of the Euler 
technique. With GEARIB, time steps varied from 0.1 second at the beginning of the 
solution up to 225 seconds toward the end of the solution. The solution time was 
less than 10 percent of the Euler solution time and less than 25 percent of the 
Crank-Nicholson solution time. 


EFFECT OF MODELING ON ALGORITHM PERFORMANCE 

This section gives some insight concerning how modeling details can affect the 
performance of transient solution algorithms - especially explicit algorithms. The 
structure chosen for the study is an insulated cylindrical shell shown in figure 12. 
The cylinder is 720 inches in length and 180 inches in diameter. The aluminum is 
0.1 inch thick and the insulation (insul) is 2.0 inches thick. The outer surface of 
the insulation is heated over a region which consists of one-third the length and 
one-half the circumference. 

Three finite-element models are used in the study. Because of symmetry, only 
half the cylinder is modeled in each case. In model I, solid (K81) elements are used 
exclusively - 39 along the cylinder length, 4 around the circumference, and 3 through 
the depth (2 elements in the insulation and 1 in the structure). The outer surface 
has quadrilateral elements which receive the heating and radiation elements (R41) 
which radiate to space. Model I contains 800 grid points and 650 elements. The use 
of solid elements to model the metal layer is not necessary and was used here as an 
extreme example of the effect of thin metal elements on stiffness. In model II the 
solid elements in the metal layer are replaced by quadrilateral elements (K41) in 
which temperatures do not vary through the thickness? this is generally a good 
assumption for thin metal structures. Model II has an extra layer of solid elements 
in the insulation in order to preserve the number of grid points in the model at 
800. In model III the grid points are identical to those of model II, but the 
insulation is modeled with one-dimensional conductors (K21). This model neglects 
lateral heat conduction but, as mentioned previously in connection with the coarse 
Shuttle wing model, this effect is small for the class of insulated flight structures 
of primary interest in the present work. 

Another aspect of the effect of modeling is comparison of results from finite- 
element and lumped-parameter models. To investigate this, the MITAS lumped-parameter 
computer program (ref. 26) was used to analyze the cylinder. The finite-element 
model I was converted to a lumped-parameter model by use of the CINGEN program 
(ref. 27). The resulting lumped-parameter model contained 625 nodes as conpared to 
800 grid points in the finite-element model. The unknown MITAS temperatures are 
located only at the centroids of each lump. 

The first 2000 seconds of the cylinder temperature response to the heating rate 
history shown in figure 13 were computed in each model. The explicit (Euler) and 
implicit (backward-difference) algorithms were used for all models, and GEARIB was 
used for the SPAR models. Solution times are summarized in table 3. Model I is 
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extremely stiff as evidenced by the small time step (t s ) of 0.06 second required for 
stability of the explicit algorithm. The high stiffness is due to the use of K81 
elements to model the metal layer. In model II, the stiffness has been essentially 
eliminated by replacing the three-dimensional elements in the metal layer with two- 
dimensional elements. In model II, the explicit technique uses time steps ranging 
from 2.4 to 10 seconds and is faster than backward-difference and GEARIB tech- 
niques. In model III, due to low stiffness again, the explicit algorithm is faster 
than the implicit but GEARIB is slightly faster than the explicit technique. 

MITAS computation times are shown in table 3. Because none of the SPAR models 
is equivalent to the MITAS model in terms of the number of unknown temperature or 
nodal connections, no direct comparison of MITAS and SPAR solution times are 
appropriate. The MITAS model is not particularly stiff as evidenced by the large 
time step used in the explicit solution technique. SPAR models II and III which 
begin to resemble the MITAS model in certain respects are also less stiff and favor 
explicit algorithms . 

Figure 14 contains temperature histories of a point in the metal portion of the 
cylinder 45° inside the edge of the heated region and 18 inches from the end of the 
cylinder. Model II is considered to be the best of the models (recall the additional 
insulation elements used); thus, the temperatures represented by the dashed line are 
expected to be the most accurate. These results are bracketed by results from 
model I and MITAS (from above) and by model III (from below). There are negligible 
differences between temperatures from the implicit and explicit solutions for any 
given model. Results from models II and III are different from those of model I 
because of the extra layer of insulation elements. The MITAS temperature history 
agrees well with that of model I (on which the MITAS model is based) except for some 
differences beginning at about 1400 seconds. 

Figure 15 contains temperature distributions from each model along a generator 
of the cylinder 45° from the edge of the heated region which contains both the heated 
and unheated parts. Temperatures are plotted for the aluminum surface of the 
cylinder. All four models give consistent results and are able to track the large 
temperature drop across the boundary between the heated and unheated regions. Of 
particular note is the close agreement of the results from models II and III. This 
close agreement tends to verify the adequacy of using one-dimensional conduction 
elements in the insulation of model III, which led to substantial solution time 
savings. The results of this modeling study show that modeling can have a major 
impact on the stiffness of the system and on solution times, especially for an 
explicit algorithm. 


EFFECT OF MATRIX SCALING 

For models I and II of the cylinder problem (table 3), the GEARIB algorithm, 
despite using much larger time steps, was only marginally faster than the Crank- 
Nicholson implicit algorithm, and for model II it was slower than the explicit 
method. This is partially due to the different ways of handling the temperature- 
dependent material properties. In the explicit and implicit methods as implemented 
in SPAR, the properties are represented as piecewise constant within time intervals 
specified by the user (by the input quantity TI) in SPAR. Material properties for 
each element are evaluated at the beginning of each interval based on the average 
temperature of the element nodes, and the conductivity and capacitance matrices are 
regenerated only at those times. Results for models I, II, and III in table 3 were 
obtained by using TI = 20 seconds. For the GEARIB algorithm, the material prop- 
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erties are assumed to vary continuously and the residual R (eq. (7)) must be 
evaluated at each iteration used in solving equation (8). Thus, there is a need for 
continuous regeneration of the conductivity and capacitance matrices. This extra 
effort is the price paid for higher accuracy. Furthermore, this burden of recal- 
culating matrices and residuals shows up most noticeably in problems which utilize 
solid (K81) elements because of the size of the matrices for those elements 
(model III does not contain K81 elements). A way to ease the burden (for thermally 
isotropic elements) has been identified and has been recently implemented in SPAR. 
The method is to generate the matrices only once for unit values of the appropriate 
property and simply scale the matrices by the property whenever it is updated. The 
effect of scaling on the solution time is shown in table 4. The most significant 
impact of scaling is in the implicit and GEARIB solution times in models I and II 
which contain K81 elements and the least impact is seen in model III which has 
none. Scaling thus helps to ease the burden on implicit algorithms of regeneration 
of conductivity matrices for the K81 elements and makes them more attractive for use 
in large complex models. Scaling has a minimal effect on the explicit solution time 
for model I because most of that time is due to time marching and relatively little 
is due to matrix regeneration. 


PARTITIONING METHODS 

The examples in the previous section have shown that stiffness of the problem is 
the major factor in choosing an explicit or an implicit temporal integration method 
for thermal analysis. When the problem is stiff, an implicit method is preferable; 
when the problem is nonstiff, an explicit method is more efficient. Some problems, 
however, call for a mix of implicit and explicit techniques. These problems occur 
when the high stiffness is due to only a part of the model. This typically happens 
when a small number of elements with high diffusivity or short conduction lengths are 
used. Recently, several promising techniques have been developed which combine the 
attributes of explicit and implicit integration techniques. These techniques are 
referred to as partitioning techniques and are still in the research stage in 
contrast to the more mature status of the purely explicit and implicit integration 
techniques. A survey of partitioning methods may be found in reference 28. This 
section presents examples of the use of such methods for transient thermal analysis. 


The Hughes-Liu Mixed Explicit-Implicit Technique 

One class of partitioning methods is based on dividing the model into explicit 
and implicit parts. The Hughes-Liu mixed implicit-explicit technique, which was 
originally developed for problems with fluid structure interaction (refs. 17 and 18), 
is an example of such a partitioning method. The basic idea is to use implicit inte- 
gration for the portion of the model which is responsible for the stiffness of the 
problem and explicit integration for the rest of the model. This idea was adapted by 
Malkus, Reichmann, and Haftka (ref* 19) to transient heat transfer in a structure 
which is composed of a good insulator and a good conductor. The structure is divided 
into an explicit group of elements and an implicit group of elements. Equation (1) 
is then rewritten as 


T - T 

2C “ + Kl(T n + T n +1 ) + ~ Q(T n } " S(T n+1 ) = 0 (9) 

n 
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where K is the partition of the conductivity matrix corresponding to the elements 
in the explicit group, and K 1 is the partition for the elements in the implicit 
group. When all the elements are explicit, the method reduces to the explicit Euler 
method; when all the elements are implicit, the method reduces to the implicit Crank- 
Nicholson method. The method was applied to the configuration in figure 16 which is 
a coarse model of one-half the configuration in figure 3. The insulation finite 
elements are in the explicit group; the metal finite elements are in the implicit 
group. The diffusion ratio between the insulation and the metal was kept constant at 
50, whereas the absolute values of the diffusivity were changed to control the stiff- 
ness of the problem. Table 5 contains two sets of results comparing the calculation 
times for the fully explicit (Euler) , fully implicit (Crank-Nicholson) , and Hughes- 
Liu techniques for two values of the dif fusivities . For the lower value of the dif- 
fusivity, the problem is not stiff and the explicit method performs better than the 
implicit method. In this case, the Hughes-Liu method is comparable in performance to 
the explicit method. For the higher value of the diffusivity (higher stiffness), the 
implicit method is faster than the explicit method because of the stiffness -*f the 
problem. In this case, the Hughes-Liu technique is significantly better than either 
the explicit or the implicit method. This simple example indicates the potential of 
the Hughes-Liu technique for improving the efficiency of transient thermal analysis 
for more complicated insulated- structure problems. 


Operator-Splitting Algorithms for Linear Equations 

The Hughes-Liu technique represents partitioning at the level of integration of 
the differential equations. Other partitioning methods may be appropriate later in 
the solution process. An important class of these methods are used in the solution 
of the linear algebraic equations generated by Newton's method utilized in implicit 
solution techniques (eqs. (7) and (8)). For notational simplicity, the equations 
generated by Newton ' s method are written as 


Ax = b (10) 


Operator-splitting methods are based on separating the matrix A into two parts M 
and N such that 


A = M - N 


(ID 


Equations (10) and (11) are combined and the resulting equation is solved 
iteratively. A simple fixed- point iteration is used 


x m+1 = M -1 (b + Nx 111 ) 


( 12 ) 


This iteration may be modified to include an overrelaxation parameter to and is 
written as follows : 


x m+ 1 = COM - 1 (b + Nx 1 ") + (1 - cojx 1 " 


(13) 
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The choice of M and N is based on the following general considerations. First, 
the matrix M should be significantly easier to invert than the matrix A. Second, 
it should be chosen for fast convergence of the iteration in equation (13). The 
choice of M may also be based on matrix topology considerations - for example, the 
sparsity of the matrices A and M - or it may be based on the properties of the 
mathematical model that led to matrix A. Consider the latter choice first. 

Operator splitting based on model characteristics .- If the problem is composed 
of a good conductor and a good insulator, the splitting (which parallels the Hughes- 
Liu technique) is to include in M the contribution of the good conductor and the 
diagonal term from the good insulator. This choice was applied to the analysis of a 
simple metal insulation system (shown in fig. 17) by David S. Malkus, Elwood T. 

Olsen, Peter I. Reichmann, and Raphael T. Haftka under NASA Grant NSG-1266 at 
Illinois Institute of Technology. This model depicts an insulated slab heated by 
internal sources and cooled at the outer insulation surface by radiation to free 
space. The Crank-Nicholson algorithm with the modified Newton's method was used to 
generate the equations to be solved. Computation times with and without operator 
splitting were compared for several grids. The results are given in table 6. Though 
the splitting was effective (from the standpoint that the solution of eq. (13) 
required only a single iteration) , the results in table 6 are somewhat disappoint- 
ing. The results indicate this type of operator splitting is effective only if 
frequent factoring of the Jacobian is needed. 

Operator splitting based on matrix structure .- Several choices for the 
matrices M and N based on the matrix topology are available. Two obvious choices 
that result in an easily invertible M matrix are the diagonal of matrix A and the 
lower triangular part of A. These choices lead to the well-known Jacobi and Gauss- 
Seidel iterative methods, respectively. The Gauss-Seidel iteration with over- 
relaxation is commonly known as successive overrelaxation ( SOR) . Although matrices 
based on the above are easily inverted, they usually result in slow convergence of 
the iteration process. When the matrix A is very sparse and has a large bandwidth 
(as when it is generated by a three-dimensional finite-element model) , a choice of 
M which is based on a partial elimination or incomplete Cholesky decomposition that 
retains the sparsity of A has been suggested. In this case, 

M = LL T (14) 


where L is a lower triangular matrix which has the same sparsity pattern as A. 

This incomplete Cholesky decomposition has been combined with the conjugate gradient 
method by Meijerink and Van der Vorst (ref. 20) to produce the incomplete Cholesky 
conjugate gradient ( ICCG) algorithm. The conjugate gradient method forms the basis 
for other iterative algorithms , such as preconditioned conjugate gradient methods 
(e.g. , ref. 29). The algorithm has been applied by Haftka and Kadivar (ref. 30) to 
the insulated cylinder problem. For this problem four methods were compared: 

(1) Gaussian elimination (a direct solution technique); (2) SOR; (3) incomplete 
Cholesky decomposition with fixed-point iteration (eq. (12)); and (4) the ICCG 
algorithm. Results are shown in figure 18 for cylinder models having 400, 720, and 
1100 grid points. These results demonstrate the advantages of iterative methods over 
elimination, especially for poorly banded problems. Also, the ICCG algorithm is 
consistently superior, for the problem studied, to the two other iterative procedures 
and is nearly independent of the matrix bandwidth but depends on the degree of 
sparsity. 
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CONCLUDING REMARKS 


This paper discussed an effort to obtain increased efficiency in calculating 
transient temperature fields in complex aerospace vehicle structures. Explicit 
solution techniques which require minimal computation per time step and implicit 
techniques which permit larger time steps because of better stability were 
reviewed. A promising set of implicit solution algorithms having variable time steps 
and order, known as GEARIB, was described. Test problems for evaluating the 
algorithms were defined and finite- element models of each one were described. The 
problems included a coarse model of the Space Shuttle Orbiter wing, an insulated 
frame test article, a metallic panel for a thermal protection system, and detailed 
models of single-bay and three-bay sections of the Shuttle wing. Calculations were 
carried out using the SPAR finite- element program. Results generally indicated that 
implicit algorithms, particularly the GEARIB techniques, are more efficient than 
explicit algorithms for solution of transient structural heat-transfer problems when 
the governing equations are stiff. Stiff equations were encountered in most of the 
test problems and are frequently encountered in thermal analysis of insulated thin 
metal structures . 

Studies were made of the effect on algorithm performance of different models of 
an insulated cylinder test problem. These studies revealed that the stiffness of the 
problem is highly sensitive to modeling details and that careful modeling can reduce 
the stiffness of the resulting equations to the extent that explicit methods are 
quite effective. 

Evaluations of two partitioning techniques were also performed. First, a mixed 
implicit- explicit technique adapted to thermal analysis was demonstrated. In this 
method, the model is separated into stiff and nonstiff portions with the explicit 
algorithm applied to the nonstiff part and the implicit algorithm applied to the 
stiff part. This method was comparable in performance to the explicit method for a 
nonstiff problem and faster than either the explicit or implicit techniques for stiff 
problems. Second, two operator-splitting techniques for speeding up the solution 
of the algebraic equations associated with implicit algorithms were described and 
demonstrated. Both are based on separating the coefficient matrices into two parts 
and solving the resulting equations by iteration. The first technique bases the 
separation on stiff and nonstiff elements of the structure, and the second bases 
the separation on the sparsity structure of the matrix. The most effective of the 
techniques tried was a separation based on matrix sparsity (using an incomplete 
Cholesky decomposition) with conjugate gradient iteration. This technique was 
especially noteworthy in that its performance appeared to be insensitive to the 
band structure of the matrix. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
July 1, 1982 
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TABLE 1.- PERFORMANCE OF ALGORITHMS FOR TRANSIENT THERMAL ANALYSIS 

OF VARIOUS MODELS 

[Solution time is CPU time] 


Explicit algorithm 

Implicit algorithm 

Euler 

Crank-Nicholson 

GEARIB 

Time 

step, 

sec 

Solution 

time, 

sec 

Time 

step, 

sec 

Solution 

time, 

sec 

Time 

step, 

sec 

Solution 

time, 

sec 

(a) Space Shuttle Orbiter wing (45 00- sec temperature history) 

10 

2 288 

10 

11 730 

1 to 528 

570 

(b) Shuttle frame (2000-sec temperature history) 

0.16 

1 723 

1 

10 

25 

50 

475 

249 

106 

65 

50 to 170 

54 

(c) Multiwall thermal protection system (3200-sec temperature history) 

0.007 

a 98 400 

1 

5 

b 28 400 
6 352 

1.0 to 113 

2754 

(d) Three-bay section of Shuttle wing (3500- sec temperature history) 

2.6 

10 560 

5 

16 800 

0.1 to 229 

1950 

(e) Single-bay section of Shuttle wing (3500-sec temperature history) 

0.1 

3 205 

1 

1 145 

0.1 to 225 

245 


a Extrapolated value based on 12 296 seconds for 400 seconds of 
temperature history. 

^Extrapolated value based on 8879 seconds for 1000 seconds of 
temperature history. 
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TABLE 2.- EFFECT OF TIME STEP ON ACCURACY OF IMPLICIT 
ALGORITHMS FOR SHUTTLE FRAME 


Time 

step, 

sec 

Temperature of 
node 309 a at 
1200 sec, °F 

Temperature of 
node 49 a at 
1200 sec, °F 

1.0 

335.7 

398.6 

10.0 

335.6 

398.5 

25.0 

331.6 

396.0 

50.0 

328.3 

394.7 

b . 16 

335.7 

398.6 

c 50 to 170 

337.5 

400.3 


a See figure 3(b). 
^Explicit algorithm. 
C GEARIB . 


TABLE 3.- EFFECT OF MODELING ON SOLUTION TIME FOR INSULATED 

CYLINDER PROBLEM 

[Solution time is CPU time] 



Explicit algorithm 

Implicit algorithm 




Crank-Nicholson 




Euler 

and backward- 

GEARIB 

Model 



difference 




Time 

Solution 

Time 

Solution 

Time 

Solution 


step. 

time , 

step. 

time. 

step. 

time. 


sec 

sec 

sec 

sec 

sec 

sec 

SPAR model I 

0.06 

10 107 

10 

1880 

1.0 to 83 

MM 

SPAR model II 

2.4 to 10 

1 518 

10 

1920 

5 to 106 

WE&m 

SPAR model III 

3.3 to 10 

279 

10 

536 

2 to 133 

■PM 

MITAS lumped- 






WUmM 

parameter model 

10 

226 

10 

320 


MM 
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TABLE 4.- EFFECT OF MATRIX SCALING ON SOLUTION TIME FOR 
INSULATED CYLINDER PROBLEM 

[Solution time is CPU time] 



Explicit algorithm 

Implicit algorithm 


Euler 

Crank-Nicholson 

GEARIB 

Model 














Time 

Solution 

Time 

Solution 

Time 

Solution 


step. 

time. 

step. 

time. 

step. 

time. 


sec 

sec 

sec 

sec 

sec 

sec 

SPAR model I 







With scaling 

0.06 

8 000 

10 

735 

1 to 83 

520 

Without scaling 

0.06 

10 107 

10 

1880 

1 to 83 

1779 

SPAR model II 







With scaling 

2.4 to 10 

392 

10 

782 

5 to 106 

505 

Without scaling 

2.4 to 10 

1 518 

10 

1920 

5 to 106 

1707 

SPAR model III 







With scaling 

3.3 to 10 

247 

10 

504 

2 to 133 

231 

Without scaling 

3.3 to 10 

279 

10 

536 

2 to 133 

266 
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TABLE 6.- EFFECT OF OPERATOR SPLITTING ON SOLUTION TIME FOR INSULATED SLAB 
[Solution time is CPU time on Prime 400 computer] 


Mesh 

Number of 
time 
steps 

Number of 
factorings 

Solution time 

, sec, for - 

Modified Newton's 
method without 
splitting 

Modified Newton's 
operator method 
with splitting 

X 

X 

91 

20 

119 

149 

7x7x1 

52 

26 

95 

99 

7x7x1 

51 

51 

137 

117 

X 

— i 

X 

t> 

91 

20 

237 

340 

7x11 x 1 

51 

51 

234 

223 

7x7x2 

91 

20 

1205 

1129 

7x7x2 

52 

26 

1164 

752 

7x7x2 

51 

51 

1911 

928 














(a) Configuration. 



(b) Finite- element model. 

Figure 1.— Space Shuttle Orbiter wing. 
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Figure 2.- Transient temperatures in Shuttle Orbiter wing 
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Figure 5.- Multiwall thermal protection system panel 
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Figure 6.- Transient temperatures in titanium multiwall TPS panel. 



Bay 1 Bay 2 Bay 3 
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Figure 8.- Finite-element model of three-bay section of Shuttle wing (wing station 328). (From ref. 25. 




Figure 9.- Heating rates on model of three-bay section of Shuttle wing. (From 











Heating rates on single bay of Shuttle wing. 
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.Figure 12.— Finite-element models of insulated cylinder. 
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in parentheses are values of DT in seconds 
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(a) Configuration and boundary conditions. 



(b) Finite-element model. 


Figure 16.- L-shaped insulated metal configuration used to evaluate 
Hughes-Liu mixed implicit-explicit integration algorithm. 


Internal heat 



Figure 17.- Insulated slab problem used to assess operator splitting. 
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Computation time, CPU sec 



0 10 20 30 40 50 60 70 80 

Half-bandwidth 

(a) 400 grid points. 


Figure 18.- Comparison of various linear operator- splitting techniques in 
transient thermal analysis of insulated cylinder. (CPU time is for 
Prime 400 computer. ) 
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Figure 18.- Concluded. 
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